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ABSTRACT 

Relativistic Sweet-Parker type magnetic reconnection is investigated by relativistic resistive magne- 
tohydrodynamic (RRMHD) simulations. As an initial setting, we assume anti-parallel magnetic fields 
and a spatially uniform resistivity. A perturbation imposed on the magnetic fields triggers magnetic 
reconnection around a current sheet, and the plasma inflows into the reconnection region. The inflows 
are then heated due to ohmic dissipation in the diffusion region, and finally become relativistically hot 
outflows. The outflows are not accelerated to ultra-relativistic speeds (i.e., Lorentz factor ~ 1), even 
when the magnetic energy dominates the thermal and rest mass energies in the inflow region. Most of 
the magnetic energy in the inflow region is converted into the thermal energy of the outflow during the 
reconnection process. The energy conversion from magnetic to thermal energy in the diffusion region 
results in an increase in the plasma inertia. This prevents the outflows from being accelerated to 
ultra-relativistic speeds. We find that the reconnection rate TZ obeys the scaling relation TZ ~ 5"°'^, 
where S is the Lundquist number. This feature is the sam e as that of non- re lativistic reconnection. 
Our results are consistent with the theoretical predictions of lLvubarskvl ()2005f ) for Sweet-Parker type 
magnetic reconnection. 

Subject headings: magnetic fields - magnetic reconnection - magnetohydrodynamics (MHD) - rela- 
tivistic processes 



1. INTRODUCTION 



Magnetic reconnection is one of the most impor- 
tant subjects in the st udies of space, laboratory, and 
astro physical plasmas (jBiskampI 119861 : iPriest fc Forbes! 
120001 ). In particular, it plays an essential role in the 



understanding of energy conversion processes in high 
energy plasmas that characterize astrophysical com- 
pact objects, such as neu tron stars (Kennel & C oronitil 
119841: iLvubarskv &: Kir k 2001; Lvubarskv 2003) soft 
gamma-ray repeater s (iLYutikoy 2006; Gill fe Hevl<l2010l; 
Masada et al.l 120101 ) . active galactic nuclei (iDi Matted 
19981). and gamma- ray b ursts ( iDrenkhahnI 120021: 



McKinnev fc Uzdenskvl[2010t iZhang fc Yan 20Tj 



While issues concerning the physical mechanisms 
and properties of magnetic reconnection remain un- 
settled in the relativistic regime, a few theoreti- 
cal studies of re l ativist i c effects have been made 
(iBlackman fc Fieldl IT991 iLvutikov &: Uzdenskvl [200l 



Lvubarskv! 1200 5': 'Tolstvkh et al.' 120071: iTenbarge et alj 
20101 ). !Lvu'tik ov & Uzdensky (2003) studied relativis- 



tic Sweet-Parker type reconnection within the frame- 
work of magnetohydrodynamics (MHD), finding that the 
reconnection-driven outflow can have an ultra-relativistic 

|takahashiQcfca.jp| 



speed (Lorentz factor 7^1) when the magnetic energy 
is preferentially converted to kinetic energy. They con- 
cluded that the reconnection rate would be enhanced in 
the relativistic regim e due to Lorentz contraction. In 
contrast, [Lyubarski^ (|2005! ) concluded that the outflow 
cannot be accelerated to a relativistic speed (7 ~ 1) be- 
cause the magnetic energy should be converted into ther- 
mal energy to maintain the pressure balance across the 
current sheet. The effect of the Lorentz contraction is 
then negligible and the reconnection rate would not be 
enhanced. 

Recently, further numerical studies have been con- 
ducted on relativistic magnetic reconnection. Particle-in- 
cell simulations are mainly employed to ascertain the re- 
connection mechanism in the collisionless regime without 
introducing a phenomenological parameter, i.e., electric 
resis tivity (Zenitani & Hoshino 2001: Jaroschek et aL| 
2004; Zenitani & Hoshino 2007; Zenitani & Hesse 2008i 
^essho & Bhattach ariee, 2007: Zenitani fe Hesse . 20Q8bl). 
Watanabe & Yoko vama! (|2006! ) studied Petschek type 
relativistic magnetic reconnection for the first time in 
a spatially localized resistivity model using relativistic 
resistive magnetohydrodynamic (RRMHD) simulations. 
Re lativistic two-fl u id MH D si mulations were perfor med 
by !Zenitani et all ()2009a! fbl) . !Zenitani et aD pOlOl ) nu- 
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TABLE 1 
List of Simulation Runs 



model 


/9o 


O"0 


"A,0 


Rai 


Tmax 




t20 (1 - 




""out 






7out 


B1R2 


0.1 


20 


2.00 


200 


5.0 


2.45 


191.4 


22.7 


1.86 X 10" 


-2 


33.7 


1.40 


B2R2 


0.2 


10 


1.41 


200 


3.0 


2.00 


195.1 


13.8 


6.58 X 10" 


-3 


18.7 


1.29 


B4R2 


0.4 


5 


1.00 


200 


2.0 


1.73 


206.5 


9.36 


2.22 X 10' 


-3 


11.1 


1.17 


B8R2 


0.8 


2.5 


0.71 


200 


1.5 


1.58 


233.3 


7.15 


8.27 X 10" 


-4 


7.48 


1.08 


B2R1 


0.2 


10 


1.41 


100 


3.0 


2.00 


290.7 


13.4 


9.56 X 10- 


-3 


19.5 


1.20 


B2R5 


0.2 


10 


1.41 


500 


3.0 


2.00 


157.6 


14.1 


2.21 X 10" 


-3 


13.9 


1.28 


B2R10 


0.2 


10 


1.41 


1000 


3.0 


2.00 


151.9 


14.4 


1.00 X 10" 


-3 


9.91 


1.22 



Note. — Columns: (6) possible maximum outflow Lorenz factor evaluated from equation |6]l; (7) Ma, max = 
i^max/^^A.o ; (§) timc at which the current sheet length reaches Y — 20 = L20] (9) ^l^^ total specific enthalpy 
evaluated at (X, Y) — (5, 0) when t — £20; (10)-(12) the magnetization parameter, the specific enthalpy, and the 
Lorcntz factor evaluated at (X, Y) — (0, L20) when t — t2Q. 



merically examined relativistic magnetic reconnection 
using various resistivity models and obtained Petschek 
and Sweet-Parker type magnetic reconnection in rela- 
tivistic plasmas. 

While these previous studies have revealed much about 
Petschek type reconnection in the regime of low mag- 
netic Reynolds number, Rm ^ 160, the physics of Sweet- 
Parker type magnetic reconnection and the dependence 
of the reconnection rate on the magnetic Reynolds num- 
ber remain unsolved. In this Letter, we focus on Sweet- 
Parker type magnetic reconnection and investigate the 
basic properties through RRMHD simulations. This is 
the first systematic study of relativistic Sweet-Parker 
type magnetic reconnection for a broad range of mag- 
netic Reynolds number Rm < 10^. 

2. NUMERICAL MODEL 

We numerically solve a set of RRMHD equations in 
two-dimensional Cartesian coordinates {X,Y) with the 
simple form of Ohm's law. 



J =qv 



V X B — {E ■ v)v] 



(1) 



where q, j, 77, v, 7, E, B are the charge density, electric 
current, electric resistivity, three velocity, Lorentz factor, 
electric field, and magnetic field, respectively. We set the 
light speed, Boltzmann constant, and average particle 
mass as unity throughout this paper. 

The r elativistic Harris sheet is adopted as an initial 
setting (jKirk fc Skiaeraasenll2003[ ). The spatial distribu- 
tion of the magnetic field is then given by Bx{X,Y) = 
B^{X,Y) = and By{X,Y) = tanh(2X/A), where 
Bq is the field strength of the sheath plasma {X — >■ ±00) 
and A is the initial thickness of the current sheet. The 
gas pressure and density profiles are determined by the 
local pressure balance, i.e. p = pQ + [BQ—By{X, y)]/(87r) 
and p = p/Tq, where po is the initial gas pressure of the 
sheath plasma and Tq is the initial plasma temperature, 
which is assumed to be constant throughout the region. 
We set po = To = 1 and fix the specific heat ratio T as 
F = 4/3 throughout this paper. The magnetic reconnec- 
tion is triggered around the origin by a perturbation of 
the magnetic field described in vector potential form as 
5A^ = -5Bo\cyi^[-{X^ + Y^)/\\ where 5Bo = O.OSSq 
is the amplitude of the perturbation. 

The plasma /3 in the sheath is defined as /3o = 87rpo / . 
The magnetization parameter corresponding to each 
model is tJo = 20, 10, 5, 2.5 , where a = B^ / {An p-/^) 
and a subscript "0" denotes a quantity of the sheath 



plasma. The corresponding Alfvcn four speed is ua,o = 

VAfll — v\ Q = 2.0, 1.41, 1.0, 0.71, where va is the 

Alfven speed. In our models, the electric resistivity is 
assumed to be uniform. We vary the resistivity rj for 
magnetic Reynolds number Rm = AttX/t] — 100, 200, 
500, 1000 . We note that Rm defined by the full thick- 
ness of the initial current she et is twice as large as that 
used in the previous study bv lZenitani et al.l ()2010[) who 
defined Rm by the half thickness of the current sheet. 
The parameters adopted in each model arc summarized 
in Tabic □ 

The two-dimensional calculation is performed in the 
X-Y plane with a volume bounded by X = [0, 50] and 

Y = [0, 75]. The length and time are normalized by the 
initial thickness of the current sheet A and its Alfven 
crossing time ta = X/va,o- We use a non-uniform grid 
in the X-direction and a uniform grid in the y-direction 
of 790 X 3000 zones. The minimum grid size in each 
direction is Ax = 5 x 10~^ and Ay = 2.5 x 10~^. A 
symmetric boundary condition is applied at X = and 

Y = 0. A free boundary condition is imposed at X = 
50 and Y = 75. Wc^calculate numerical fluxes using 
the HLL method (|Harten et ai.lll983l : [Komissarovll2007f l 
with an operator-splitting method. We use the implicit 
scheme to solve Ampere's equation in order to maintain 
numerical stabilit y when the electric resistivity is small 
(jPalenzuela et all [2009) . The implicit scheme enables 
us to study magnetic reconnection with larger magnetic 
Reynolds number than previous studies. 

Before presenting the numerical results, we estimate 
the speed of outflow driven by magnetic reconnection. 
Consider steady Sweet-Parker type magnetic reconnec- 
tion in a diffusion region of length 26 in the inflow direc- 
tion and 2L in the outflow direction. Mass and energy 
conservation between the inflow and outflow give 



Pin'"in2i — PoutUout2'5, 



(2) 



Pm 



T -r 



Pin 7inUin + 



Pout 



r- 1 



An 

Pout I 7outWout H 



2L 

An 



26, (3) 



where u = is the four velocity. The subscripts "in" 
and "out" indicate physical quantities in the inflow and 
outflow regions. Note that the ideal MHD {E — —v x B) 
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Fig. 1. — Snapshots at i = 230 of models B2R2 (Figs, la, lb, and Ic) and B2R10 (Fig. Id). The colors show (a) the outflow component 
of the four velocity {uy), (b) the plasma temperature {T = p/p), (c) the electric current perpendicular to the X—Y plane (jz), and (d) the 
gas pressure (p). The solid lines show the magnetic field lines. The arrows in Fig. la represent the velocity fields. 



can describe both inflow and outflow plasmas because 
these flows are outside the difl'usion region. Combining 
equations (2) and (3), we obtain 



(1 + (Tin + /lin) 7in = (1 + ^out + Kut) 7o 



(4) 



where h — rp/[(r — l)p] is the specific enthalpy. This 
is similar to the Bernoulli equation that describes the 
conservation of the total enthalpy flux between the inflow 
and outflow. Assuming Uout ^ 1, which is reasonable for 
the case of anti-parallel {B^ = 0) magnetic reconnection, 
we obtain 

l + ^in+fin fr\ 

Tout = —T 7in. (5) 

1 + /lout 

If the magnetic energy is converted preferentially into 
kinetic energy (h-m = feput), as was ideally assumed in 
iLvutikov fc Uzdenskvl ()2003l ). we obtain an upper limit 
of the outflow Lorentz factor 7i„ax, 



7n 



1 



-7ii 



(6) 



listed in Table [1] is evaluated from an initial a-m 
and h (we take 7in = 1 in all models for numerical eval- 
uation). Equation ^ reduces to 7max = (1 + o'in)7in 
when the thermal energy is negli gible (see also equation 
32 in ILvutikov fc Uzdenskvl[2003| ) . The outflow becomes 
super- Alfvenic in this case. 

We stress that 7max is the upper limit of the outflow 
Lorentz factor under the condition hin — hout- Since a 
part of the magnetic energy should be spent for plasma 
heating by ohmic dissipation, the enthalpy of the outflow 
is expected to increase and its Lorentz factor would be 
less than 7max except in the ideal situation. 

3. RESULTS 

Figure [1] shows snapshots at t — 230 of model B2R2 
(Figs, la, lb, and Ic) and B2R10 (Fig. Id). We fo- 
cus here on Fig. la-lc (model B2R2), and will refer 
to Fig. Id in § m The colors show (a) the outflow 
component of the four velocity {uy), (b) plasma tem- 
perature (T = p/p), (c) electric current perpendicular 
to the X-Y plane (jz), and (d) gas pressure (p). The 



solid lines depict the magnetic field lines and the ar- 
rows represent the velocity fields. Following the initial 
perturbation, the magnetic field lines start to reconnect 
around the origin. As the reconnection proceeds, the 
current sheet is elongated along the ±y-direction, indi- 
cating the formation of a Sweet-Parker current sheet. 
The plasma temperature inside the current sheet Tcs is 
almost uniform with a constant value (~ 8.7Tin). If adi- 
abatic heating is the main process for plasma heating, 
Tcs would increase up to ~ {pcs/ PinY^^Ti^ ~ 1.22Tin, 
where p^s is the density inside the current sheet. Here 
we used the numerical result that pout — Pes ^ l-Spin- 
Since the evaluated temperature is much less than the 
observed value, the plasma should be heated by ohmic 
dissipation rather than adiabatic heating in the diffusion 
region. The heated plasma is accelerated in the diffusion 
region along the ±y-direction, resulting in the forma- 
tion of hot outflows. The maximum outflow speed at 
this time is ^ 0.66, which is slightly smaller than the 
Alfven speed of the sheath plasma {va,o — 0.82, dis- 
cussed later). The reconnection outflows collide with a 
magnetic bubble (plasmoid) that originates in the initial 
Harris sheet. This causes the formation of reverse fast 
shocks around y ~ 35 outside the diffusion region. The 
plasma temperature near the plasmoid increases because 
of shock heating and the plasma expands. We note that 
a backflow (vy < 0) structure forms around the plasmoid. 
Such backflow structures are also observed in Pctschek 
type magnetic reconnection, w hich would be cause d by 
the expansion of the plasmoid ([Zenitani et al.|[2010[ ). 

At this stage, the electric field around the reconnec- 
tion region is well developed. The typical strength of the 
electric field near the reconnection point {X ^ Y ^ Q) 
is Ez ~ 1.0 X 10~^i?o. In addition, an electric field 
{Ez — VyBx) is also induced near the plasmoid by the 
reconnected magnetic fields that are swept by the out- 
flow and accumulate around Y = 37. Its maximum am- 
plitude is Ez — 0.17-Bo) which is 17 times larger than 
that near the reconnection point. Such strong electric 
fields near pl asmoids are also observed in particle-in-cell 
simulations (jJaroschek et al.l 120041 : iZenitani fc Hoshinol 
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Fig. 2. — Time evolution of the maximum Alfven four Mach 
number along X = (Ma)- The diamonds denote Ma at t = 420- 



I2007D and two-fluid simulations (jZenitani et al.ll2009bl ). 
These electric fields can generate high-energy particles 
(IHoshinol l2005r). 

Figure [2] shows the time evolution of the maximum 
outflow component of the four velocity normalized by 
the Alfven four speed in the sheath plasma along X = 0, 

Ma = max[uy{0,Y)]/uA,o- Each line shows the time 
evolution of a different model listed in Table [T] After the 
onset of magnetic reconnection, the outflow is accelerated 
along the F-direction. Its acceleration decreases as the 
reconnection develops (e.g., t > 200 for model B2R2). 
We note that Ma does not increase to the upper limit 
of the Alfven Mach number Ma, max = ""max/uA.o listed 
in Table [T] for each model, where Umax = \/7max ~ 1- 
This suggests that the magnetic energy is not efficiently 
converted into kinetic energy. The saturation values of 
Ma are independent of Rm, but weakly decrease as ctq 
increases. 

In order to estimate the energy composition of in- 
flow and outflow plasmas, we evaluate the magnetization 
parameter a, specific enthalpy h and total specific en- 
thalpy 1 + a + h at {X,Y) = (5,0) for inflows and at 
{X, Y) = (0, 20) for outflows, summarized in Table [TJ 
These values are evaluated when the length of the cur- 
rent sheet L, which is deflned by jz(0, L)/jz(0, 0) = a, 
reaches L — L20 = 20 (hereafter, we refer to this time as 
^20)- We take a — 0.5 but the results are independent 
of a when a < 1 since jz is almost constant inside the 
current sheet and rapidly decreases at the edge of the 
current sheet. For example, the length of the current 
sheet of model B2R2 at t = 230 is L = 33.3 (see Fig. 

Since we take ctq > 1, the inflow plasma is dominated 
by the magnetic energy (1,/iin ^ o-m)- In the outflow 
plasma, /lout is much larger than CTout and unity, i.e., 
the thermal energy density dominates the rest mass and 
magnetic energy densities. Moreover, /lout is comparable 
to the total specific enthalpy of inflow (1 -I- (Tin + /iin)- 
These results show that almost all the magnetic energy 
is converted into thermal energy due to joule heating. 
By using this fact, we can evaluate the outflow Lorentz 
factor from equation ([5]) as 7out — 1 in the condition 
7in — 1. The outflow is not accelerated to a relativistic 
speed 7max- This result is consistent with that of the an- 
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Fig. 3. — Time evolution of the reconnection rate (7?.). The 
diamonds denote 7J. at t = 420- 



alytical work bv lLvubarskvl ()2005D . They concluded that 
the magnetic energy should be converted into thermal en- 
ergy ((Tin — /iout) to maintain the pressure balance across 
the current sheet. We confirmed that the current sheet 
reaches the pressure balance in the reconnection system, 
resulting in the formation of relativistically hot outflows. 
Such hot outflows cannot be accelerated to a relativistic 
speed 7max since the magnetic energy is spent in heating 
rather than acceleration. Moreover, a larger enthalpy 
{h > 1) leads to an increase in the plasma inertia in 
the relativistic regime (see equation (3). This also pre- 
vents the outflow from being accelerated. Sub-Alfvenic 
outflow is a natural outcome of relativis tic Sweet-Parker 
reconnection because 7out — \/ Uout + 1 remains of order 
unity when ua,o increases with (Tq. Thus, the Alfven four 
Mach number decreases with increasing ctq in the rela- 
tivistic regime. This conclusion is different from that for 
non-relativistic reconnection, in which the outflow speed 
reaches the Alfven speed. 

Figure [3] shows the time evolution of the reconnection 
rate TZ — —Vx/va at {X,Y) = (5,0). An initial rapid 
increase in the reconnection rate is caused by the ini- 
tially imposed perturbation. After that, the reconnection 
rate gradually increases with time and reaches a maxi- 
mum. Although the outflow speed is almost maintained 
at the saturated value, the reconnection rate decreases 
after passing its peak because the current sheet contin- 
ues to elongate. Such an elongation of th e current sheet 
is als o observed in non-relativistic plasma (jTanuma et al.l 
l2001f ). According to the non-relativistic theory, the as- 
pect ratio of the Sweet-Parker current sheet is described 
hy S/L = TZ. Since the curvature radius of the initial 
magnetic field is infinite in the Harris sheet (i.e., uniform 
in the ^-direction), L tends to increase as the reconnec- 
tion proceeds, resulting in a decrease in the reconnection 
rate. 

Figure 2] shows the reconnection rate TZ at {X, Y) — 
(5, 0) as a function of Rm at t = ^20- The reconnection 
rate is almost independent of ctq for a fixed Rm = 200, 
while it decreases with Rm for a fixed tro = 10. The 
solid line in this figure s hows the relation TZ = S""^/^, 
obtained by iLvubarskvl ()2005D (see equation 8 in that 
paper), where S = 4:TtL2ova,o/v = Rm{L2o/ ^)yA,o is the 
Lundquist number. In this plot, we assumed <tq = 10 to 
evaluate the Lundquist number that weakly depends on 
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Fig. 4. — Reconnection rate (TZ) as a function of magnetic 
Reynolds number (Rm)- The reconnection rates are evaluated 
at {X,Y) = (5,0) when t = t20. The black diamonds corre- 
spond to models of different ijjv/ with the same (tq = 10, while 
the blue square, green asterisk, and red cross correspond to models 
of (To = 20, 5, 2.5 with the same Rm = 200, respectively. The 
solid line shows TZ = S~^'^, where S is the Lundquist number, 
proportional to Rm- 



(To- 

We found that the reconnection rate is well fitted by 
the steady model TZ = S~^-^, while the system is not 
exactly steady state. Although we are not sure why our 
time-dependent results follow those of the steady model, 
we also found the relation ti/trcc = 0.5 holds in the 
Sweet-Parker regime (e.g., 170 ^ t < 240 for the model 
B2R2), where tree = —{d\ogTZ/dt)^^ is the decreasing 
time of the reconnection rate and tj^ = {d\ogL/dt)~^ is 
the increasing time of the length of the current sheet. 
This relation is obtained from TZ = S'~"-^ by assuming 
that TZ and L are time dependent. This means that the 
reconnection system evolves while maintaining the rela- 
ti on TZ = S-°-^. 

iLvutikov fc Uzdenskvl ()2003[ ) proposed that the recon- 
nection rate is enhanced by the relativistic effect of the 
Lorentz contraction by a factor when ao ^ 1 (see 
equation 39 in that paper). In our numerical models, 
however, such effects never become effective because of 
mildly relativistic outflow (7 ~ 1), as discussed above. 
Therefore, the reconnection rate is almost independent 
of (To and is described by S~°-^. We conclude that rela- 
tivistic Sweet-Parker reconnection is a slow process for 
energy conversion, as for non-relativistic plasma. 

4. CONCLUSION AND DISCUSSION 

We have developed a relativistic resistive magnetohy- 
drodynamic (RRMHD) code that is applicable to plas- 
mas with larger magnetic Reynolds numb ers than pos- 
sible in previous stu dies (jWatanabe fc Yok ovama 2006: 
iZenitani et al.l 120101 ). We confirmed that the reconnec- 
tion outflow does not accelerate to a relativistic speed 



because the magnetic energy released by magnetic re- 
connection is spent on plasma heating rather than accel- 
eration. The plasma heating results in increasing iner- 
tia, which also prevents the outflow from being acceler- 
ated. Since the Lorentz factor of the outflow is of order 
unity, the enhancement of the reconnection rate due to 
the Lorentz contraction is ineffective. Thus, the recon- 
nection rate of the relativistic plasma obeys the relation 
TZ = S^'^'^, which is the same as that for non-relativistic 
plasmas. We confirmed that this relation holds in the 
Sweet-Parker regime (i.e., tio < t < t^o, where tio and 
^30 are the time at which L = 10 and L — 30, respec- 
tively). Thes e results are consi stent with the theoretical 
prediction of iLvubarskvl (|2005[ ). 

We note that humps appear in the reconnection 
rate beyond the Sweet-Parker regime (see e.g., t > 
200 for the model of B2R10 in Fig. E]). At this 
stage, we observed a growth of the tearing instabil- 
ity in the elongated current sheet (Fig. Id, the 
model of B2R10). An increase in the reconnec- 
tion rate following the instability is also observed 

in non-relativistic 

20011: iBhattachariee et al.l 
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reconnection 
[2?M 



Jc 



Shibata & Tanurna 
Cassak et al., 200!; 



120101 : iHuang fc Bhattacharieel i2010j . 
Once the tearing instability develops, the current sheet is 
disrupted and the system evolves to a non-steady state. 
A remarkable feature is that the increase in the outflow 
velocity coincides with the enhancement of the recon- 
nection rate (e.g., t > 200 of model B2R10 in Figures [2] 
and|3]). Although the reason for the acceleration is not 
yet clear, we speculate that the outflow is accelerated 
by the pressure gradient force. When the plasmoid cre- 
ated due to the tearing instability flows along the cur- 
rent sheet, the plasma in the current sheet is swept up 
by it. The subsequent outflow can then be accelerated 
by the pressure gradient force (see Fig. [TJl). Since the 
plasma inertia h ^ p/p decreases while it is accelerated 
by the pressure gradient force, the outflow speed might 
be ultra-relativistic j ^ 1. In such cases, relativistic 
effects would facilitate the energy conversion in the re- 
connection process. These scenarios need to be verified 
through numerical simulations and will be reported in a 
subsequent paper. 
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